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ABSTRACT 

We compute the rates of capture of stars by supermassive black holes, using time-dependent Fokker-Planck 
equation with initial conditions that have a deficit of stars on low-angular-momentum orbits. One class of initial 
conditions has a gap in phase space created by a binary black hole, and the other has a globally tangentially- 
anisotropic velocity distribution. We find that for galactic nuclei that are younger than ~ 0.1 relaxation times, 
the flux of stars into the black hole is suppressed with respect to the steady-state value. This effect may 
substantially reduce the number of observable tidal disruption flares in galaxies with black hole masses M, > 
10^ M0. 

Subject headings: galaxies: nuclei — galaxies: kinematics and dynamics 


1. INTRODUCTION 

A Star passing at a small enough distance from a super- 
massive black hole (SMBH) is captured or tidally disr upted, 

S cing a detectable flare in multiple wavebands (iReesI 
. Estimating the rate of such events is the topic of 
the loss-cone theory, first developed in 1970s in applica¬ 
tion to hypothetic al intermediate-mass black holes in glob¬ 
ular c lusters (e.g. iFrank & Reeslll9'^ iLightman & Shapir^ 
IT97^ . Later on, this theory was applied to SMBH in 
galac t ic centres dS ver & U lmenll99^lMagorrian & Tremainel 
Il999t IWang & Merritli l2004li . typically estimating the tidal 
disruption rate to be in the range 10“"* — 10“® events per 
year per galaxy. Observ ational constraints from optical 
ivan Velzen & Far^l2014ft. UV (iG ezari et akl l2008h o r X- 

rav (iDonlev et al.l 120021 lIGiabibullin & SazonovI 120141) sur- 

veys roughly fa ll into the same range (see a l so a r eview by 
lKomossdl2015h . Recently IStone & Metzger! (1201 4ft raised a 
concern that the observationally derived event rates are sys¬ 
tematically lower than the theoretical estimates, by as much as 
one order of magnitude, although the uncertainties are hardly 
smaller than that on either side. 

In this paper we investigate a mechanism that can signif¬ 
icantly reduce the event rates in galaxies with long enough 
relaxation times and alleviate the tension between theory and 
observations, namely the influence of tangentially anisotropic 
initial conditions. The source of this anisotropy could be a 
gap in the low-angular-mom entum region of the pha se space, 
created by a binary SMBH (iMerritt & Wangl 12005ft . or sim¬ 
ply a mild bias towards circular orbits. We do not attempt to 
model individual galaxies, but focus instead on the compari¬ 
son of steady-state rates, typically used in the literature, with 
the reduced ones from anisotropic but otherwise the same ini¬ 
tial conditions, in order to derive the suppression factor as a 
function of galaxy parameters. We only consider spherically 
symmetric galaxies, therefore obtaining a lower boundary on 
the possible event rates. 

The paper is organized as follows. In Section |2] we write 
down the Fokker-Planck equation that describes the diffusion 
of stars in angular momentum, and derive its time-dependent 
analytical solution. Then in Section |3] we present the modi- 
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fications of initial conditions that have tangential anisotropy. 
In Section |4] we obtain the estimates of capture rates for our 
choice of initial conditions, and compute the suppression fac¬ 
tor with respect to the steady-state capture rates. We discuss 
the implications of our results in Section |5] 

2. ANALYTICAL SOLUTION OF THE DIFFUSION EQUATION 

Two-body relaxation that changes energies E and angu¬ 
lar momenta J of stars can be descri bed in terms o f orbit- 
averaged Fokker-Planck equation ('e.g. lMeiTTtil20133. Chap¬ 
ter 5). It is commonly assumed that the changes in an¬ 
gular momentum are much more important for computing 
the c apture rates (iFrank & ReesIIl976t ILightman & Shanira 
Il977ft . and for this reason the diffusion in energy is usu¬ 
ally neglected; numerical solution of the more general two- 
dimensional F okker-Planck equation (iCohn & Kulsrudiri978l: 
iMerritll 12015ft confirms the validity of one-dimensional ap¬ 
proximation for processes that occur on a timescale much 
shorter than the relaxation time Uei- At a fixed energy E, 
the distribution function (DF) of stars /(j) as a function of 
normalized angular momentum j = J/ JdrdE), where Jdrc 
is the angular momentum of a circular orbit with the same en¬ 
ergy, satisfies the diffusion equation in a cylindrical geometry: 

dfiE,j,t) Vdfdf\ 

dt Aj dj y dj J ■ 

Here D(£'~) oc tC.} is the orbit-averaged diffusion coefficient 
(e.g. lMeri^l2013bl equation 18), which is assumed to be in- 
depe ndent of ?, allowing an analyt ical solution to this equa¬ 
tion (iMilosavlievic & MerrittlUOO^ . The boundary condition 
at J = 1 is of the Neumann type: df /dj — 0 (zero-flux 
condition). The presence of the black hole creates a cap¬ 
ture boundary at j = jic, the angular momentum at which 
a star would be tidally disrupt ed at periapsis. As discussed 
in ILightman & Shapiro (11977ft . there are two limiting cases 
for the behaviour of f{j) near j\c, depending on the ratio 
q = VT/j/^ between the mean-square change in j due to 
relaxation over one orbital period T and the size of the loss 
cone. In the general case, one may use a Robin-type bound¬ 
ary condition (a linear combination of the function and its 


















































2 


derivative), which naturally interpolates between the regimes 
of empty {q <C 1) and full {q 3> 1) loss cone: 


/(jlc) 


«jic df{j) 

2 dj 


3 — 3\c 


a{q) 


(g" + g4)i/4. (2) 


In the steady state, the flux of stars towards the capture 
boundary T = {'I2/A)j df{j)/dj is nearly independent of 
j at low j, and the solution of Equation [T] has a logarithmic 
profile: /(j) c>c logj + const. Extrapolating this functional 
form to a ll j, one obtains the class ical quasi-steady-state solu- 
tion (e.g.lVasiliey & Merrittl2013l equation 51, with TZ = p). 
ICohn & Kulsrudid 19781) defined an effective capture boundary 
jo = jic exp(—a/2), at which the inward extrapolation of the 
logarithmic solution reaches zero; this gives the same bound¬ 
ary condition at the true capture boundary as Equation |2] ex¬ 
cept for their slightly different expression for a. 

The general time-dependent solution of equation ([T]l subject 
to the speci fied boundary conditions can be expressed in a 
series form (iMilosavljevic & Merrittil200^ equations 24-26): 


fU,t) = CmA{l3m,j) exp(-I?/3^f/4) , (3) 

m—1 

A{p,j) = uP3)ym - YomJiW), ( 4 ) 


c„ 


JUPtuAc) - JKPm) 


j'A{l3m,j') f{j',0)dj\ 


where Ji and Ij are Bessel functions of the first and the sec¬ 
ond kind, A are the basis function, /3m are the roots of a cer¬ 
tain equation that satisfy the boundary conditions for each 
basis function, and Cm are the expansion coefficients com¬ 
puted from the initial conditions /(j, 0). The above study 
derived their expressions for the limiting case of an empty 
loss cone, i.e. the boundary condition /(jic) = 0. Tak¬ 
ing t he same inward extrapo lation of the logarithmic profile 
as in ICohn & KulsrudI (119781) . one can use their expressions 
with a modified capture boundary jo for the general case of a 
non-empty loss cone (this method was adopte d for the time- 
dependent solution in IVasiliev & MenTt3l2013l) : however, the 
validity of this approach depends on the assumption of the 
solution being close to the steady-state profile. Instead, in 
this work we generalize the analytical solution to the arbitrary 
boundary conditions given by Equation|2] Namely, the coef¬ 
ficients /3m that enter the expressions for the series expansion 
are the roots of the following equation: 


[<70 (/3m jic) T 0:/3mjlc'7l (/3mjlc)/2] Ti(/3m) 

[YoiPmJlc) + a/3mjlcij(/3mJlc)/2] Jl(/3m) = 0. (5) 


This modification is necessary to obtain a rigorous solution 
for non-trivial initial conditions, such as those considered in 
the next section, although it gives only a minor correction 
to the more approximate treatment in the case of a nearly- 
logarithmic form of the solution (i.e. late enough into the evo¬ 
lution). 

The series solution with a finite number of terms rn-max 
breaks down at small times {t < so we com¬ 

pute the solution numerically with a finite-difference scheme 
in this case. 


3. TANGENTIALLY ANISOTROPIC INITIAL CONDITIONS 

The stationary solution with a logarithmic j-dependence of 
DE is attained at times greater than an appreciable fraction of 





Fig. I.— Monte Carlo simulations of a binary SMBH in a galaxy with 
7 = 1 Dehnen density profile. The binary mass is 10“^ of the total stellar 
mass, the mass ratio g = 1, and the binary started on a nearly-circular orbit 
at separation 0.2, roughly equal to the radius of influence. 

Top panel: Phase space (squared angular momentum vs. energy) after the 
binary has cleared the low angular momentum region. Dashed green and 
d ot-dashed blue lines sho w the definition of the gap region in this study and 
in IMerritt & Wand 420050 : solid red line marks the angular momentum of a 
circular orbit. 

Bottom panel: approximation of the initial distribution function (Equation[6) 
used in this work. 


the relaxation time. On the other hand, this time may well ex¬ 
ceed the age of the Universe, especially in low-density galac¬ 
tic nuclei. Therefore, the choice of initial conditions becomes 
important for determining the present-day capture rates. We 
consider two classes of initial DEs with a deficit of stars at 
low angular momenta, compared to the isotropic population. 

The first possible reason for such a deficit is the ejection 
of stars by a binary SMBH that may have existed in a galaxy 
previously. The slingshot mechanism is responsible for the 
flattening of the density profil e and the form ation of cores 
(IMilosavljevic & Merritll 120011: iMeiritll l2006h. which have 
been detected ob servationally (e.g. iDullo & Grahamll201^ 
[Rush et ^l2013h . More importantly, it creates a gap in the 
phase space, ejecting stars with angular momenta less than 
some critical value, which corresponds to the periapsis radius 
comparable to the radius of a hard binary, «h = <?/[4(l + 
7 )^] Tinfl, where q = m 2 jmi < 1 is the mass ratio of the bi¬ 
nary, and Tinfl is the SMBH radius of influence. In this work 
we adopt the definition of rinfl as the radius containing the 
mass of stars equal to 2 (toi -I- m 2 ) before the slingshot pro- 
c ess has started. 

iMerritt & Wangl (120051) considered a time-dependent so¬ 
lution for the one-dime nsional Eokker-Planck equation , us¬ 
ing the expressions of IMilosavljevic & Merritll (120031) . i.e. 
an empty loss cone boundary condition, and taking the ini¬ 
tial distribution in angular momentum as a Heaviside step 
function: f{E,J) (x 0(J — Jgap(i3)). They defined the 

gap width as Jgap(E) = Kah\/2, [E — fl)(ifah)], with a 
dimensionless constant iT ~ 1. We have used a modi- 
fled version of the Monte Carlo code Raga (IVasilievll2015l: 
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IVasiliev et al.ll2015h to determine the distribution of stars in 
angular momentum in a galaxy with a binary SMBH; our 
results are better fit by an energy-independent gap width, 
Jgap = a/ K' G(mi -I- TO 2 )ah, with K' ~ 3, and a more 
gradual drop towards smaller J (Figure[T]i: 


f{E, J, 0) = f{E) ■ min j j . ( 6 ) 

Another possible choice of initial conditions involves a DF 
that has a tangential anisotropy at all J, not just a gap at 
small J. The simplest possibility is to consider DF in a 
factorized form; f{E,J) = {I - I3)[J/ f{E), 

where /(£’) is the counterpart of the usual isotropic DF and 
is computed f r om a given density profile with the method 
of ICuddefordI (119911) . Such DF corresponds to a constant 
velocity anisotropy coefficient /3 (iBinnev & Tremain3l2008l 
equation 4.61); the case of weak tangential anisotropy /3 = 
— 1/2 is consistent wit h some observationa lly-based models 
of galactic centers (e.g. lThomas et al.ll20l4 figure 2), and re¬ 
sults in a very simple expression for DF: 


f{E) 


JcircjE) (f r p($) 

Stt^ r($) 


(7) 


4. RESULTS 

We have considered a set of iDehnenI (119931) 7 -models with 
a black hole mass M, = 10“^ of the total mass in stars. The 
energy-dependent part of DF f{E) or f{E) and the diffusion 
coefficient T>{E) are computed numerically from the given 
density profile, using the Eddington inversion formula or its 
Cuddeford’s generalization. We explore several values of the 
power-law index 7 of the central density profile, and for each 
value of 7 we chose to consider a one-parameter family of 
models by scaling M, and rjnfl simultane ously, according to 
the following relation (iMerritt et al.ll200^ : 

ri„fl = ro[M./10®MQ]0'^6. (8) 

As our default normalization, we set tq = 30 pc, but we also 
consider values of tq — 20 and 45 pc. The ratio of the in¬ 
fluence radius to the scale radius of Dehnen profile is {0.091, 
0.047, 0.016} for 7 = {0.5,1,1.5}. 

It is natural to express our results in dimensionless 
units: the flux normalized to the steady-state capture 
rate, and the time meas ured in units of relaxation time 
(IBinnev & Tremainel2008l equation 7.106) at nnfl. For 10® < 
M,/Mq < 10® and 0.5 < 7 < 1.5, these scale approxi¬ 
mately as 

Ig [E^t/ (Mgyr"^)] « -4.6 - 1.5 lg(ro/30 pc) (9) 
+ 0.2(1- 7 ) lg(M./lO®M 0 ), 


Ig [frei/yr] - 13 -f 0.4(1 - 7 ) + 1.5 lg(ro/30 pc) (10) 
-f 1.28 lg(M./10®MQ). 

The steady-state flux is comparable to M./frei- Figure |2] 
top panel, shows the ratio of time-dependent to steady-state 
capture rate, as a function of time, for two representative mod¬ 
els with a gap. We also plot the same quantity measured at the 
time 10 ^° years, for various choices of 7 , M, and g; in this 
manner, the abscissa corresponds to the black hole mass ac¬ 
cording to our scaling (fTOl i. 



10 “'’ 10"^ 10^ 10^ 10° 


i/^rel 


Fig. 2.— Suppression factor — the ratio of time-dependent to steady-state 
capture rate, as a function of time normalized to the relaxation time at the 
radius of influence. Individual points correspond to models with different 
density profiles and M* = 10 ®, 10 ® ®,..., 10 ® Mq; the abscissae corre¬ 
spond to the Hubble time (10^® yr) measured in units of relaxation time, thus 
the most massive black holes are at the left side of the plot. 

Top panel: families of models with an initial gap in angular momentum dis¬ 
tribution due to the action of a pre-existing binary SMBH; different symbols 
encode 7 ; open symbols are for models with smaller gap width (binary mass 
ratio q = 1 / 10 ) and filled — to equal-mass binaries; half-filled symbols are 
for q = 1 and different normalizations of rinfi—M* relation. Dotted and dot- 
dashed curves show examples of time-dependent flux for particular choices 
of parameters. 

Bottom panel: tangentially anisotropic (/3 = —1/2) models for different 7 . 


The time time required to establish the stea dy-state profile 
at a given energy is ~ ( Jgap/Jcirc)^ irei (e.g. lMerrit l l2013bt 
equation 34); as the maximum of the total flux arrives from 
energies corresponding to rinfl, the time to refill the gap is 

roughly frefill ~ (ah/Anfl) frel(nnfl), Or 


frefill ^ 10 ^® yr X 


q 

4(l-f g)2 


M, rp 

10^ MqJ V30pc 


1.5 

( 11 ) 


For t < freflii the flux is reduced compared to the stationary 
value, which is commonly used in calculations of tidal dis¬ 
ruption rates. The maximum value of capture rate reached at 
t ^ frefiii is somewhat lower than the steady-state value, due 
to the fact that the J-averaged DF is also depleted at high 
binding energies (where Jgap > Jcirc) with respect to the 
value used in the steady-state calculation. The flux reaches 
1/2 of its maximum value at ti /2 — O.lfrefiii- Moreover, at 
t > trp fiii it starts to decline in the absense of diffusion in 
energy. iMerritli (120151) has performed numerical integration 
of two-dimensional (E, J) Fokker-Planck equation restricted 
to the region inside rjnfl, also using initial conditions with a 
gap at J < Jgap, and found a qualitatively similar behaviour 
if the diffusion in energy was artificially switched off. On 
the other hand, taking it into account modifies the solution at 
t > 0.1/rei SO that it tends to a steady-state profile. Therefore 
we may trust our calculations roughly up to a time when the 
flux reaches its maximum. 
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We also explore the effect of changing the normalization 
in — M, relation (l8]l. This, of course, modifie s both the 
time-dependent flux and the relaxation time at rinfl (1911 Oi l, but 
the normalized values still stay on the same curve for each 
7 and q. Finally, the second class of models with globally 
tangentially anisotropic initial conditions (Figure |2] bottom 
panel) produce a milder decline of the capture rate at f < 
O.lfrel. 

5. DISCUSSION AND CONCLUSIONS 

We have considered the question whether the capture rate 
of stars by SMBH can be substantially lowered with respect 
to the steady-state value by a suitable modification of initial 
conditions. We used an analytical time-dependent solution of 
the Fokker-Planck equation to compute the capture rate for 
the given initial conditions as a function of time. Two classes 
of initial conditions were analyzed: a gap in the low-angular- 
momentum region of phase space, created by a binary SMBH, 
and a separable distribution function with a mild tangential 
velocity anisotropy. 

The results of our study can be summarized as follows. If 
the relaxation time in the galaxy centre is short enough (< 
10 ^^ yr), any difference between initial conditions is erased 
quite rapidly, and the capture rate approaches the steady-state 
value. On the other hand, for tre\ ^ 10^^ yr anisotropic ini¬ 
tial conditions may substantially reduce the capture rate (in 
the intermediate range of frei, their effect is moderate). In our 
series of models with a phase-space gap, for SMBH masses 
M, > 10 ^ Mq the suppression factor (the ratio of time- 
dependent to steady-state capture rates) drops quite rapidly, 
plunging below 10“^ for M, > 10® Mq. Note that for even 
heavier black holes, visi ble flares constitute a small fraction 
of the captured stars (e.g. lMacLeod et aTll2012L figure 15). In 
models with mild tangential anisotropy the suppression factor 
is not so extreme, but nevertheless can reduce the flux by a 
factor of a few for M, > 10® M(^ _ 

We have extended the work of iMerritt & Wand (120051) into 
the range of smaller SMBH masses (< 10® Mq), since they 
are more numerous in the Universe and are expected to dom¬ 
inate the overall tidal disruption rates. Moreover, smaller 
M, residing in more compact galactic nuclei are not well 
described by the empty-loss-cone regime adopted in that pa¬ 
per. In this study we have derived the solution for the general 
case; we checked that assuming an empty-loss-cone bound¬ 
ary condition does not substantially change the results for 
M, > 10^ Mq, but overestimates the flux by a factor of a few 
for the least massive SMBHs. We used a somewhat different 
initial distribution function inside the gap than the above pa¬ 
per, but checked that adopting their initial conditions changes 
the results only marginally. However, our estimates of the 
time required for the capture rate to reach 1/2 of its steady- 
state value, ti /2 — O.lfrefiib the latter given by Equation [Til 
are about an order of magnitude longer than shown in fig¬ 
ure 3 or equation 14 of iMerritt & Wan j (l2005l) for the same 
galaxy parameter^ We believe that this discrepancy might 
be due to a calibration error in th at paper, as our express ions 
for /refill agree with equation 7 in iMerritt & Szelll (l2006h and 


equation 36 in iMerritj (l2013bl) . 

We deliberately have made a number of simplify¬ 
ing assumptions that drive our capture rates towards 
lower values. First, we assumed a spherical ge¬ 
ometry, while it is known that non-spherical torques 
can r pult in a higher efflciency of loss-cone repopula¬ 
tion ([Ma gorrian & Tremain j [19991: iMerritt & PoonI 120041 : 
Holley^ ockelmann & Sigurdssonll2006l) . Naturally, the dif¬ 
ference between spherical and non-spherical geometry starts 
to manifest itself above the same threshold value of a s 
the influence of initial conditions (figure 4 in lVasilie^l2014l) . 
again underlining the distinction between relaxed and non- 
relaxed galactic nuclei. Note that by relaxed we here mean 
the systems that had enough time to establish a nearly steady- 
state logarithmic profile in angular momentum distribution; 
this does not mean they were able to relax in energy space as 
well and develop a lBahcall & WolfI (119761) cusp. Second, we 
neglected non-cla ssical phenomen a that may increase the re- 
laxation rate (e.g.|Alexande3l2012l) . such as mass segregatio n 
(iFreitag et al.ll200^ massi ve perturbers (|Perets et alJl20d^ . 
or resonant relaxation (e.g. lMerrittjl2013aL Chapter 5.6): the 
latter, however, is typically not very effective in boosting 
the capture rates (iHopman & Alexandeil 120061). A number 
of oth er refinements have been shown bv IStone & Metzg^ 
(120141) to have little impact on the final values. Binary 
SMBH th at have not yet c oalesced also suppress the cap¬ 
ture rates (IChen et al.ll2008l) . although they may demonstrate 
brief epis odes of increased encounter rates at early stage of 
eyolu tion (llyanoy et al.l2005UChen et al.l201 IblWegg & Bodel 
1201 ll) . Thus our findings can be regarded as robust lower lim¬ 
its on the capture rates for single SMBH. 

The effect inyestigated in this paper is unlikely to substan¬ 
tially reduce the rate of tidal disruption eyents for black hole 
masses smaller than ^ lO’^ Mq. Although the yolumetric rate 
of tidal disruption eye nts is dominated by galax ies with small¬ 
est M, (e.g. figure 8 in iStone & Metzgeii2014ft . the rate of ob¬ 
servable eyents (figure 10 in that paper) may be significantly 
affected by the suggested mechanism if either (1) the SMBH 
mass function is suppressed at low-mass end, or (2) flare emis¬ 
sion mechanisms are inefficient for low-mass SMBH. These 
factors are still the main sources o f uncertainty in the esti - 
mates of the obseryable eyent rate (IStone & Metz get! 120141) : 
howeyer, it is notable that the distribution of actually obseryed 
flares peaks around M, = 10^ Mq (figure 11 in that paper). 
Thus the effect considered here may potentially be important 
in comparing the theoretical predictions of tidal disruption 
rates with obseryations. 
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